Numerical computations of facetted pattern formation in snow crystal growth 
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Facetted growth of snow crystals leads to a rich diversity of forms, and exhibits a remarkable sixfold symmetry. Snow crystal 
structures result from diffusion limited crystal growth in the presence of anisotropic surface energy and anisotropic attachment 
kinetics. It is by now well understood that the morphological stability of ice crystals strongly depends on supersaturation, crystal 
size and temperature. Until very recently it was very difficult to perform numerical simulations of this highly anisotropic crystal 
growth. In particular, obtaining facet growth in combination with dendritic branching is a challenging task. We present numerical 
simulations of snow crystal growth in two and three space dimensions using a new computational method recently introduced by 
the authors. We present both qualitative and quantitative computations. In particular, a linear relationship between tip velocity 
and supersaturation is observed. The computations also suggest that surface energy effects, although small, have a larger effect 
on crystal growth than previously expected. We compute solid plates, solid prisms, hollow columns, needles, dendrites, capped 
columns and scrolls on plates. Although all these forms appear in nature, most of these forms are computed here for the first time 
in numerical simulations for a continuum model. 
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1. Introduction 

Snow crystals grown from a supersaturated vapour lead to a 
variety of complex and often very symmetric patterns. Crys- 
tallisation from vapour is a fundamental phase transition, and a 
good understanding is crucial for many applications. Numerous 
experiments have been performed, and compilations of pho- 
tographs of artificial and natural snowflakes reveal their beauty 
and complexity, see 1 1 ] and the review [2]. The precise forms 
of snow crystals depend in a very subtle way on the temperature 
and the supersaturation. Nakaya Q analysed these dependen- 
cies in detail and combined his observations in his now famous 
Nakaya snow crystal morphology diagram, see Figure [T] At 
temperatures just below the freezing temperature thick plates 
grow at lower supersaturations and plate-like dendritic forms 
appear at higher supersaturations. At temperatures around 
-5°C solid prisms grow at lower supersaturations and hollow 
columns and needle-like crystals at higher supersaturations. If 
the temperature is decreased below -10°C one observes thin 
solid plates at low supersaturations, whereas dendrites form at 
high supersaturations. Below -25°C again columns form at 
high supersaturations. The results from Nakaya [lj, which led 
to the snow crystal morphology diagram, have been confirmed 
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by many subsequent experimental studies. Although the exper- 
iments give a rather clear picture, the physics behind the snow 
crystal morphology diagram are not yet understood. 

A continuum mathematical modelling of snow crystal growth 
leads to a quasi- static diffusion problem for the diffusion of 
the vapour molecules. The diffusion equation has to be solved 
together with rather complex boundary conditions on the free 
boundary between vapour and solid. The conditions on this 
interface are given by the continuity equation relating the flux 
of vapour molecules onto the interface to the interface veloc- 
ity, and an equation describing the attachment kinetics taking 
surface energy effects into account. In the latter condition the 
hexagonal anisotropy of snow crystals also enters. Due to being 
highly nonlinear, and since it is geometrically very involved, 
the complete free boundary problem is difficult to analyse the- 
oretically. However, there exists a large literature on numerical 
computations for diffusion limited growth and the formation of 
dendrites, which we now briefly discuss. 

Numerical approaches for crystal growth usually employ ei- 
ther sharp interface models, in which the solidification front is 
tracked explicitly, or phase field models, in which the solidi- 
fication front is modelled by a thin diffusional layer. In sharp 
interface approaches the front is described with the help of a 
parameterisation, see (3l[4l[5][6]l, or by using a level set func- 
tion, see (7). In a phase field method a new order parameter 
-the phase field- is introduced, which at the interface rapidly 
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Figure 1 : The Nakaya diagram illustrates which snow crystal forms appear at 
different temperatures and supersaturations. This figure is taken from \2\. 



changes its value between two fixed values, which describe the 
different phases, see (8j [9l [TCJ [TT]|. A popular discrete model 
for the simulation of crystal growth are cellular automata, see 
[13J H31. Moreover, molecular dynamics simulations are 
used to understand the surface structure of ice 1 15 ]. Although 
many computations have been performed, a quantitative numer- 
ical description of facet growth in combination with dendritic 
branching is missing. 

In recent research by the authors a new parametric approach 
for interface motion has been developed |[T6l [T7). The method 
has the feature that the mesh quality of the interface approxi- 
mation, which is given by a polyhedral surface mesh, remains 
good during the evolution - most earlier approaches had to deal 
with mesh degeneracies, e.g. by re-meshing the interface ap- 
proximation. In addition, the present authors were able to in- 
clude anisotropy effects into curvature driven hypersurface evo- 
lution in a numerically stable way. This allows the method to 
compute in situations in which the anisotropy is facetted, see 



It is the goal of this paper to demonstrate that with the nu- 
merical method introduced in (16] [T71 \19\ 13 it is possible to 
compute snow crystal growth in a qualitatively and quantita- 
tively satisfactory way. We will present computations showing 
a significant number of different types of snow crystals, such as 
solid plates, solid prisms, hollow columns, needles, dendrites, 
capped columns and scrolls on plates. 

In snow crystal growth models several parameters are not 
known to a sufficient precision. In particular, the surface en- 
ergy density as a function of orientation is not known in detail. 
In addition, the condensation coefficient, which embodies the 
attachment kinetics of how water molecules are incorporated 
into the ice lattice, is not known in detail. For example, it is not 
known how the condensation coefficient depends on the crystal 
orientation, see Q and the references therein for details. 

Our numerical computations seem to suggest that the 
anisotropy in the surface energy density might have a more im- 
portant impact on snow crystal growth morphologies than pre- 
viously expected. Of course, numerical simulations alone can- 



not decide whether this is in fact the case, but a comparison 
of numerical computations with experiments might help to un- 
derstand this issue better. In particular, it might be possible to 
obtain more precise estimates for the size of the condensation 
coefficient as a function of orientation. 

We also present some quantitative results, which first of all 
show what the relatives sizes of the quantities entering the sur- 
face attachment kinetics are. In addition, we show that the tip 
velocity for a dendrite growing into a supersaturated vapour de- 
pends linearly on the supersaturation. This linear relationship 
has been already observed in experiments for growing needles, 
see l2T1i . and our computations might help to relate parameters 
in the theoretical model to experiments. 

2. A continuum model for snow crystal growth 

We consider a continuum model for snow crystal growth, 
which consists of an ice crystal growing from water vapour, 
as discussed in e.g. El [221, and non-dimensionalize it. Let c 
denote the water vapour number density in gas. The diffusion 
equation in the gas phase, see El Eq. (2)], is then 



c t -DAc = in Q + (7), 



(1) 



where £l+(t) is the domain occupied by the gas phase and D is 
the corresponding diffusion constant. The mass balance at the 
gas/solid interface T(t) gives rise to 



_ 3c , 

£>^ = feoiid 
ov 



c)*V on T(0, 



(2) 



where c so iid « 3 x 10 28 m~ 3 is the number density for ice. In 
addition, v is the unit normal to Y(t) pointing into £l+(t) and "V 
is the velocity of T(t) in the direction v. In El Eq. (3)] the term 
c*V is neglected since c <?c c so iid- Furthermore, taking surface 
tension effects and attachment kinetics into account, we require, 
compare El Eq. (23)], 



c = c s 



t (l - Sk 7 + 



J3(v) v ki 



-) on T(t). 



Here Vkm is the kinetic velocity, S = y/(c so iid K T) « 1 nm = 
10" 3 yum, where y « 0.1 Jm~ 2 represents the typical order of the 
surface tension of ice, K « 1.4 x 10" 23 JK -1 is the Boltzmann 
constant, T is the temperature and c sat = c sat (T) is the equilib- 
rium number density above a flat ice surface, which is depen- 
dent on temperature. In addition, k 7 is the anisotropic mean 
curvature which incorporates the hexagonal anisotropy of the 
surface energy density. Moreover, J3 is the condensation coeffi- 
cient, which is denoted by a in El, which depends on the orien- 
tation of the crystal via the normal v. Finally, we complement 
([T]) with the boundary condition 



C — Con 



on dQ = da+(t)\T(t), 



(3) 



where Cr, 



<?sat + ^super describes the water vapour number 
density far away from the interface. Here, for convenience, we 
choose a domain QcR rf ,^ = 2,3, with £1+ (t) c Q, that is large 



2 



enough so that boundary effects can be neglected. Moreover, 
Csuper is related to the supersaturation £ sup er by 



super — ITIH2C) ^super ? 



(4) 



with niH 2 o ~ 3 x 10" 23 g denoting the mass of a water molecule. 
We recall that the supersaturation £ S uper appears on the vertical 
axis in Figure [T] 

It remains to introduce the anisotropic mean curvature k 7 . 
Instead of a constant surface energy density, we choose y to be 
dependent on the orientation of the interface. The effect of the 
underlying crystal structure is encoded into the surface energy 
by allowing y = y(v), where as stated above, v is the unit normal 
to the solid boundary T(t) pointing into the vapour region £l+(t). 
The total surface energy of an interface T is now given by the 
surface integral 

Jfytf) ds. 

It is convenient to extend y to be a positively homogeneous 
function of degree one, i.e. we define y(p) = \p\y(p/\p\) for 
all p ^ and refer to y as y from now on. The first variation of 
the above energy can now be computed as 



i.e. j t Jf, y(v) ds = - Jf, k 7 <V ds; where • is the tangential 
divergence on T and y f is the gradient of y, see e.g. l23ll24ll25lL 
and also lfT9ll6l. 

We now non-dimensionalize the problem. As a length scale 
we choose R, which we set to be lOOyum for snow crystal 
growth. As a time scale we choose 

~ _ ff_ Csolid 

In addition, we non-dimensionalize the concentration by intro- 
ducing 

c - c sa t _ 
u = . (5) 

c sat 

Then, in terms of the new independent variables x = x/R and 
t - t/t, we obtain (on dropping the " notation for the new vari- 
ables for ease of exposition) the equations 



c sat 
^solid 



d t u - A u - 



du 



in 

on T(0, 
on T(0, 



(6) 
(7) 



where p := (D c SSLt )/(R c soM Vkin) and a := S/R. Since c sat «: 
^soiid, we simplify the first equation to 



Au = in Q+(0- 



(8) 



We will choose y and fi of order one, and hence it will be im- 
portant to specify the order of magnitude of the quantities p and 
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supersaturation £ S uper (g/m 3 ) 


0.01 


0.02 


0.05 


0.1 


0.2 


0.3 


-1° C 


0.002 


0.005 


0.011 


0.023 


0.048 


0.069 


-2° C 


0.003 


0.005 


0.012 


0.025 


0.049 


0.074 


-5°C 


0.003 


0.006 


0.016 


0.031 


0.063 


0.094 


-10° c 


0.005 


0.010 


0.024 


0.048 


0.095 


0.143 


-15° C 


0.007 


0.015 


0.037 


0.074 


0.147 


0.221 


-30° C 


0.030 


0.060 


0.150 


0.300 


0.601 


0.901 



Table 1: Values of uqq, = c sup er/c sa t depending on T and £ S uper- 

a in Taking the values of c sa t/c so iid an d Vki n from the table 
in (2| p. 866] into account, we observe that 



^sat 



^solid ^kin 



0.71 x 10 _8 sOum)- 



independently of the temperature T . Moreover, for the time 
scale f, which depends on c sat , and hence on T, we obtain a 
range from 100 s at -1° C to 1300 s at -30° C. These time scales 
seem to be realistic when comparing with the experiments re- 
ported in |2ll26l. 

For the diffusion constant of water vapour in air we take D = 
2 x 10 7 (/mi) 2 s"\ see Libbrecht El p. 866], which is valid at a 
pressure of 1 atm. With the values of S and R mentioned further 
above we obtain 



1.42 x 10" 



a = S/R* 10" 



(9) 



If not otherwise stated, we will always choose these param- 
eters in all the snow crystal growth computations described in 
Section [4] For the boundary condition we set, on recalling ([3} 
and ([5]), 

C supgr 

u - uoci := on <9Q. (10) 



c sat 



With the help of the table in 12 p. 866] we compute several 



exemplary values for the fraction in ( 10) for different values of 
the temperature T and the supersaturation £ S uper; see Table [T] In 
effect, we appear to have reduced the two parameter variation 
of the diagram in Figure [T] to the single parameter ugo, in ( [TO] ). 
However, in our numerical simulations of snow crystal growth 
we will vary both uqq and the kinetic coefficient J3. Although 
in reality not much is known about the possible shapes and de- 
pendencies of 13, it is known that j3 strongly depends on T. Thus 
varying ft in our numerical computations may be interpreted as 
simulating different (yet unknown) temperature regimes. 

3. Numerical method and anisotropics 

For the numerical results in this paper we employ the finite 
element approximation introduced by the authors in (6| [27) in 
order to approximate solutions of ([6]), ([7]), ([5]), ( [T0| ). In the 
method a uniform time step r > is employed and the evolu- 
tion of the crystal surface is tracked with the help of parametric 
meshes T h that are independent from the bulk meshes T h on 
which the approximation u h of u is computed. The scheme uses 
an adaptive bulk mesh that has a fine mesh size hf around T h 



and a coarse mesh size h c further away from it. Here hf = 



and h c = ^ are given by two integer numbers Nf > N c , where 
we assume from now on that Q = (-H,H) d . The initial para- 
metric mesh 1^(0) consists of K® vertices, and this mesh is lo- 
cally refined, where elements become too large during the evo- 
lution. 

In order to successfully model the evolution of anisotropic in- 
terface evolution laws the authors introduced a stable discretiza- 
tion in (19| [6). We now discuss how y and ft have to be chosen 
in order to model situations with a hexagonal anisotropy. In this 
paper, we choose surface anisotropies of the form 

L 

y(P) = Yj -iP'GtpiK (ii) 



where G £ e R dxd , for I = 1 
itive definite matrices, and p 

a vector in W*. We remark that anisotropies of the form (Ml 



> L, are symmetric and pos- 
(p u ...,p d ) T e 



denotes 



admit a formulation of k 7 in ([7]) which can be discretized in 
a simple and stable way, see Il91 [6) . We will now demon- 
strate that these forms of y also allow one to model a hexag- 
onal surface energy in a simple way. To this end, let l e (p) := 

[e 2 \pf + p\{\- e 2 )f = [p\ + e 2 Zt 2 Pff for e > 0. 

Then a hexagonal anisotropy in R 2 can be modelled with the 
choice 

3 n 

y(p) = jheAp) := Yj + T } ^ ' (12) 



€=l 



where R(6) = cose) denotes a clockwise rotation 

through the angle and Go € [0, |) is a parameter that rotates 
the orientation of the anisotropy in the plane. The Wulff shape 
of (12) for 6 = 0.01 and Go = is shown in Figure [2J together 
with its polar plot V := {y(p) p : \p\ = 1}. For more details on 
Wulff shapes and polar plots we refer to | 

* we 



introduce the rotation matrices R\(G) := 



In order to define anisotropies of the form (11) in R 

sine 

cos (9 

1 

cos e sin e \ 

Ri(Q) •= ( 1 o . In this paper, we consider 

cose I 



cose 

-sin# 




and 



10 o 



3 

7(0) = 7heAP) ■= URi&0)+ 4= Z l ^(°o + { -r)?)> 
2 V3£l- 3 

(13) 

which is relevant for the simulation of snow crystal growth. Its 
Wulff shape for e - 0.01 is shown in Figure [3] together with its 
polar plot. 

We note that the Wulff shape of ( [T3] ) for e — > approaches a 
prism, where every face has the same distance from the origin. 



In other words, for ( 13 ) the surface energy densities in the basal 
and prismal directions are the same. We remark that if "Wq 



denotes the Wulff shape of ( [13] ) with 6 = 0, then the authors 
in [14 ] used the scaled Wulff shape \ "Wo as the building block 
in their cellular automata algorithm. In addition, we observe 



that the choice ( 13 ) agrees well with data reported in e.g. |29j 
p. 148], although there the ratio of basal to prismal energy is 



Figure 2: Wulff shape (left) and polar plot (right) in ] 
and 6 = 0. 



1 for (l2j with e = 0.01 




A 



Figure 3: Wulff shape (left) and polar plot (right) in R 3 for (b) with e = 0.01. 



computed as y B /y p ~ 0.92 < 1. In order to be able to model 



this situation as well, we generalise the choice ( 13 ) to 



z V3 i=1 * 

(14) 

so that now y B /y p = 7tb • 

A more generalized form of ( [T2| ) and ([13]), which also fits 
into the framework ( fTl) , is given by 



jiP) = 7hex(P) + o-\p\, 



(15) 



where cr > is a fixed parameter. For the case d = 2 we show 
the Wulff shape of ( [T5] ) for cr = 1 and cr = 5 in Figure [4] We 
note that for ( p3] ) with e = and cr > 0, both in the case d = 2 
and in the case d - 3, the corresponding Wulff shape still has 
flat parts, but is now smooth with no corners and, if d - 3, with 
no edges. Equilibrium crystal shapes with these characteristics 
can be found in certain metals l30l . and it is conjectured that 
they may be relevant for snow crystals as well l3T1l . 

As discussed in Q, and the references therein, the precise 
values of ft as a function of the normal v are not known. Hence 
one issue in our computations is to understand how different 
choices of ft influence the overall evolution. First choices for 
the anisotropy in the kinetic coefficient are J3(v) = 1 and J3 = y. 
It was discussed in [ 2 ] that the value of f3 is expected to change 




Figure 4: Wulff shape in R 2 for (T5) with cr = 1 (left) and cr = 5 (right). 
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with temperature and can vary quite drastically as a function 
of the orientation. Denoting by fi B the condensation coefficient 
of the basal directions and by /3 P the condensation coefficient in 
the prismal directions, it is for example expected that the growth 
of thin plates at T - -15°C is only possible if f3 p //3 B is large. 

In order to be able to vary the kinetic coefficient J3 signifi- 
cantly, in the case d = 3 we define for later use 



and 



Pm(p) = fyx,ttf) '.= \P\ + P\ + Pi] ^ (16) 



MP) = Mp) •= [io" 2 ' (pi + p\) + pi] 1 d7) 



with £ e N. We note that in practice there is hardly any dif- 
ference between the numerical results for a kinetic coefficient J3 
that is isotropic in the x\ - x^ -plane, such as /3fl at and /3 ta n, and 
one that is anisotropically aligned to the surface energy den- 
sity, such as e.g. J3 = J3f[ at y. Hence in all our three dimensional 
numerical simulations we always choose coefficients J3 that are 
isotropic in the x\ - X2-plane, e.g. ( [T6] ) or ( 17). 

In addition, it might be the case that the condensation coef- 
ficient J3 is considerably lower in the directions normal to the 
facets. In order to model this we choose 

fihex, [Tmax - 7hex(p)]) 

/(Tmax - Tmin) , (18) 

where we fix j3 max = 10 3 and j3 m [ n = 1, and where 

Tmax := maxy hex (p) e R> and y min := min y hex (p) e R> . 
l/*l=i l^l=i 



We note that for the 2d anisotropy ( [12] ) it holds that y max = 
7hex(e~ l6 °) and y^n = yhex(e l( %~ eo) ). For more details on the 
numerical method and the anisotropies we refer to [TH1 [HI 

a. 



4. Numerical computations 

4.1. Snow crystal growth in two dimensions 

In all computations for ([6]), (|7|, ([8]), ( [T0| in this subsection we 
will, if not otherwise stated, use the parameters ([9land choose 
the surface energy anisotropy y = yh ex defined by(T2|) with e = 
0.01 and O = ^ . The rotation in the definition of the anisotropy 
is used, so that the dominant growth directions are not exactly 
aligned with the underlying bulk meshes T h . Moreover, the 
radius of the circular initial crystal seed, T(0), is always chosen 
to be 0.05. 

First of all we study what influence the curvature and the 
velocity terms in {7]), and the supersaturation in ( 10 ) have on the 
evolution of the crystal. We choose the supersaturation UQti = 
0.004 and show the results in Figure [5] One observes that the 
hexagonal structure of the crystal forms quickly and that the 
facets become unstable and break after they reached a certain 
size - a phenomenon which is observed in experiments as well, 
see 13. 

In Figures [6}|8] we plot computations with larger supersatu- 
rations, uqq = 0.01,0.04,0.2. One clearly observes dendritic 





Figure 5: (Q = (-4,4) 2 , u dQ = 0.004, y = p = y hex ) T h (t) for t = 0, 5, . 
(left) and for t = 0, 50, . . . , 500 (right). Parameters are N f = 256, N c 
K° = 16 and t = 0.1. 
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Figure 6: (£1 = (-4,4) 2 , u da = 0.01, y = (3 = y hex ) T h (i) for t = 0, 5, 
(left), and for t = 0, 50, . . . , 200 (right). Parameters are N f = 512, Af c = 
16andr = 5xl0- 3 . 
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Figure 7: (Q = (-4, 4) 2 , u da = 0.04, y = J3 = y hex ) T h (t) for t = 0, 0.5, . . . , 5 
(left), and for t = 0, 5, . . . , 40 (right). Parameters are N f = 1024, N c = K° = 64 
andT = 2.5xl0~ 3 . 





Figure 8: (Q = (-4, 4) 2 , u dQ = 0.2, y = /3 = y hex ) Y h (t) for t = 0, 0.04, . . . , 0.4 



(left), and for t = 0, 0.4, . . . , 6.4 (right). Parameters are Nf ■ 
128andr = 2.5xl0- 4 . 



2048, N c = 



5 





Figure 11: = (-4,4) 2 , = 0.004) We take 7 = 7/^,/? = 1, N f = 256, 
iV c = 4,K° = 16 and r = 0.1 on the left and plot T h (t) for f = 0, 50, . . . , 500. 
We take 7 = y iso , J3 = y hex , N f = 512, N c = 16, K° = 16 and r = 10 -2 on the 
right and plot T h (t) for t = 0, 50, . . . , 500. 



Figure 9: (Q = (-8, 8) 2 , uqq = 0.04, 7 = yhex,fl = 1) Approximations of p^V 
(black, solid), a K™ g (blue, dashed) and a k™ x (red, dashed). Here p and a are 
as in {9}. The time interval is [0, 50]. 





Figure 10: (ft = (-4,4) 2 , u dn = 0.004) We choose 7 = y hex , p = with 
parameters Afy = 256, N c = 4, K® = 16 and t = 0.1 and plot T /z (0 for f = 
0, 50, ... , 500 on the left, and we choose a = 0, p = yhex with parameters 
N f = 2048, N c = 128, K° = 1024 and r = 10" 3 and plot T h (t) for t = 0, 3 on 
the right. 



growth, which is more enhanced at larger supersaturations. In 
addition, the evolution is much faster due to the fact that more 
water vapour molecules are available. 

It is also of interest to compare the size of the terms appearing 
in ([7]). To this end, we compare the numerical approximations 
of the terms 

ffV, a*? 8 , a K ™ ax 

where *V denotes the observed tip velocity, i.e. the velocity of 
the part of the interface furthest away from the origin, k™ 8 is 
the average of \k 7 \ on the interface and K™ ax is the maximum of 
\k 7 \. For a computation with y = yhex, P = 1 and uqq, = 0.04 
we plot these values in Figure [9] It clearly can be seen that 
the curvature contribution is larger than the velocity term. In 
our other computations we only observed for supersaturations 
around 0.2 and larger that the velocity term p'V is larger than 
the average curvature term a Ky g . 

In Figure [10| we set the velocity term to zero in the left com- 
putation, i.e. p = 0, and we set the curvature term to zero in the 
right computation, i.e. a - 0. We observe that leaving out the 
velocity term only has a very minor impact on the crystal evo- 
lution. On the other hand, leaving out the curvature term has 
a drastic effect - the front becomes very unstable. This can be 
explained as follows. Growth from supersaturated vapour is un- 
stable and while the velocity term in ([7]) without the curvature 
term can dampen the unstable modes, they are still unstable on 





Figure 12: (Q = (-4,4) 2 , u dQ = 0.004, 7 as in {15} with cr = 1 (left) and 
cr = 5 (right), J3 = y hex ) T h (t) for t = 0, 50, . . . , 500. Parameters are N f = 256, 



Af c = 4, K° = 16 and t = 0.1. 



all wavelengths. Whereas, the curvature term will stabilise the 
small wavelengths and will select a fastest growing wavelength, 
irrespective of the velocity term. 

We now study the influence of the anisotropy on the evolu- 
tion. On the left of Figure 1 1 we present an evolution with 
a hexagonal anisotropy for y and an isotropic /3. In the same 
figure on the right we take y isotropic and ft hexagonal. One 
clearly observes that anisotropy in the surface energy seems to 
be important to obtain facetted growth. 

This is underlined by the next computation, where we choose 

we 



(15) for the anisotropy 7, and let ft = yh ex - ln Figure 12 



present evolutions for cr = 1 and cr = 5. These computations 
show how important the facetted anisotropy in the Wulff shape 
is in order to obtain facetted snow crystals, recall Figure]?] 

It was suggested to the authors by Prof. Libbrecht that the 
difference of J3 as a function of v might be large with a mini- 
mum in the the directions of the facet normals. We hence took 
P = fihex,L and y = yi so in Figure [13] One observes that the 



anisotropy in the condensation coefficient f3 is large enough to 




Figure 13: (Q = (-4,4) 2 , u dQ = 0.04, 7 = 7 
0, 0.5, . . . , 5. Parameters are N f = 2048, N, 



P = Phex,L) r A (r) for t = 
256 and r = 1CT 3 . 



6 
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Figure 14: (£1 = (-4, 4) 2 , y = /3 = 1) Best linear fit for the tip velocity <V 
against the supersaturation u^q. 



lead to a six folded branching structure. But it is also clearly 
visible that the anisotropy in J3 does not lead to facetted growth. 

Next we study how the tip velocity in a growing dendritic 
crystal depends on the supersaturation. So far nothing is known 
theoretically for this dependence in the case of facetted growth, 
see E, §4], even though for simpler problems, i.e. in the ab- 
sence of facetting, a vast literature exists, see e.g. l25l and the 
references therein. As can be seen in Figure [9| the tip velocity 
after some time becomes basically time-independent. This is 
true also for our other computations with different supersatura- 
tions. We observe a linear dependence between the supersatura- 
tion and the tip velocity that the evolution eventually settles on. 
To underline this qualitative behaviour, we numerically deter- 
mine the value for the nearly constant tip velocity <V for several 



values of ugn, see Figure [14] for a plot of these velocities. We 
also show the best linear fit to this data, which is given by a 
linear function with slope ^2.35. We remark that a similar lin- 
ear relationship between velocity and supersaturation has been 
observed experimentally for needles, see |2TI . 

We end this subsection with some computations, where we 
use a time dependent choice for uqq. This models changing 
physical conditions. In particular, in the first computation we 
set 

' 0.004 t e [0, 50) u [60, 200] , 
0.08 t e [50, 60) ; 

while in the second computation we set 



(19) 



0.08 te [0, 12) U [50,56] , 
0.004 te [12,50); 



(20) 



see Figure 15 for the results. In the last computation we use a 
widely varying uqq. We set UQait) = 0.2 for t e [0,0.2), then 
Udo.it) = 0.4 till t = 0.3, then u dQ (t) = 0.08 till t = 1, then 
um(t) = 0.004 till t = 10, then u m (t) = 0.08 till t = 12, then 
uao.it) = 0.004 till t = 20, and then u da (t) = 0.08 until the end, 
see Figure [T6| for the results. 

4.2. Snow crystal simulations in three space dimensions 

Also in three space dimensions we use the physically rele- 
vant parameters introduced in Section [2] see in particular ([9]), 



Figure 15: (Q = (-4,4) 2 , y = fi = y hex ) T h (t) for t = 0, 10,..., 200 (left, 
(l% and for t = 0, 2, . . . , 12, 50, 52, 54, 56 (right, j20j). Parameters are N f = 
TU24, N r = K° = 64 and r = 10" 3 . 



and choose, if not stated otherwise, the three-dimensional vari- 



ant of yhex, see ( 13 ), with e = 0.01 and #o = f^- m a U computa- 
tions with the exception of Figure 17 the initial crystal seed was 
spherical with radius 0.05. Similarly to our computations in two 
space dimensions, we observe in our three dimensional numeri- 
cal computations that the surface energy anisotropy is important 
in order to obtain facetted growth. If we do not choose the sur- 
face energy strongly facetted, then we do not observe facetted 
growth of the crystal. 

One issue in three dimensions is to understand how the pa- 
rameter J3 leads to either horizontal flat growth or to columnar 
vertical growth, which may yield e.g. solid prisms or needles, 
respectively. First of all, we attempt to compute a self-similar 
hexagonal evolution, i.e. a crystal where the basal and prismal 
facets grow with the same velocity. This is motivated by a theo- 
retical result in [32], in which the existence of self-similar evo- 
lutions of crystals, where the Wulff shape is a cylinder, was 
shown. We choose p = a - 1, uqq = 21, y = J3 = 
as in ( 14), vary the ratio yxB = 7 B /y P and observe that, upon 
starting the evolution with T(0) being a scaled Wulff shape, for 
7tb ~ 0.95 the evolution is self-similar up to discretization er- 



1, i.e. 



rors. See Figure [T7] for a computation with yrB = 0.95. 

For the remainder of the computations we fix y TB = 
we choose y = yh ex as in ( 13 ), and use the physically relevant 
parameters in ([9]). For the first such computation we set uqq = 
0.004 and J3 = 1, see Figure [T8] We can clearly see that the 



facets of the growing crystal are aligned with the Wulff shape 
of y. We also note that facet breaking occurs both in the prismal 
and in the basal directions. In this context we refer to 1261 , 
where similar facet breaking was observed in experiments. 

It is well known that the condensation coefficient J3 varies 
strongly for different orientations, depending on the meteoro- 
logical environment. In particular, the value of can differ quite 
drastically between directions which correspond to basal facet 
normals and ones which correspond to prismal facet normals, 
see (2). We hence perform different numerical computations 
for the condensation coefficients JS^ and /3 tSL \\ defined in 
and (FT]). 



16) 



We begin with a repeat of the simulation in Figure [T8j but 
now choose as kinetic coefficient J3 = /3fl a t,2 and J3 = /3fl a t,3 ; see 



Figures 19 and 20 In comparison to the evolution in Figure 18 
one observes that the smaller condensation coefficient in basal 
directions leads to flat crystals. This is related to shapes in the 
Nakaya diagram for temperature between 0°C and -3°C and 



Figure 19: (O = (-4, 4) 3 , n 5n = 0.004, y = y hex ,/3 = /3 flat , 2 ) I*(50). Parameters 
are iV/ = 128, N c = 16, ^ = 98 and r = 10 _1 . 




Figure 16: (O = (-4,4) 2 , y = J3 = y hex ) T h {t) for t = 0, 0.2, . . . , 1 (left top), 
t = 0, 1, . . . , 12 (right top) and for r = 0, 1, . . . , 30 (bottom). Parameters are 
N f = 2048, N C = K° = 128 and r = 2.5 x 10" 4 . 





J 



Figure 17: (Q = (-8,8) 3 , y = = y™ with y T B 
0, 0.1, 0.2; and T h (0.2) within Q. Parameters are N f -- 
1538 and r= 10" 4 . 



= 0.95) T h (t) for f = 
512, N c = 32, = 



J 





Figure 18: (Q = (-4,4) 3 , w 5Q = 0.004, y = y hex ,/3 = 1) T h (50). Parameters are 
N f = 128, N c = 16, ^ = 98 and r = 10" 1 . 




Figure 20: (Q = (-8,8) 3 , u dQ = 0.004, y = y hex , p = / 
50, 100, 150, 200. Parameters are N f = 256, N c = 32, K° 



t , 3 ) T h (t) for r = 
)8 andr= 10 _1 . 



between -10°C and -22°C. 

A computation with a supersaturation uga = 0.002 and 
P - P\2lVl,\ can be seen in Figure [21] In this case the conden- 



sation coefficient is larger in the basal direction and we obtain 
a solid prism, which can be found in the Nakaya diagram at 
temperatures between -5°C and -10°C and at low supersatura- 
tions. 

At higher supersaturations uga = 0.004 we obtain for J3 = 



Aaiu the results shown in Figure 22 We also give some plots 



of the rescaled water vapour density in Figure 23 We observe 
Berg's effect (33), which states that the concentration is largest 
at the edges and decreases towards the centre of the facet. It 
is believed that facet breaking occurs when the concentration 
becomes too non-uniform on the facets (34). In Figure 22 we 



observe facet breaking for the basal and prismal directions, al- 
though the breaking predominantly occurs on the basal facets. 
Choosing the condensation coefficient even larger in the 



basal directions leads to Figure 24 We observe hollow columns 
as in the Nakaya diagram between -5°C and -10°C at low, 
but not too low, supersaturations. Increasing the condensation 
coefficient in the basal directions even further, i.e. choosing 
P - Aaii,3> l ea ds to the evolution depicted on the left of Fig- 
ure [25] On the right we also display a computation on a coarser 



grid. Both results in Figure 25 lead to needle growth, which also 
appears in the Nakaya diagram. We remark that the shape on 
the right of Figure [25] is caused by numerical noise and round- 





Figure 24: (O = (-4,4) 3 , u dQ = 0.008, y = y hex , (3 = J3 m , 2 ) T h (t) for t = 
1, 2, 5, 10, 20, 30, 40, 50; and 1^(50) within O. Parameters are N f = 128, 
N c = 16, K° = 98andr = 10 _1 . 



Figure 22: (Q = (-4,4) 3 , u dn = 0.004, y = y hex , (3 = Aaii,i) I*(0 for ? = 
1, 2, 5, 10, 20, 30, 40, 50; and r^O) within Q. Parameters are N f = 128, 
N c = 16, £j = 98andr = 10 _1 . 






Figure 26: (H = (-4,4) 3 , u da = 0.02, y = y hex , (5 = /? flat3 ) I*(f) for t = 
0.05, 0.1, 0.2, 0.3; and T h (03) within O. Parameters are N f = 512, N c = 32, 
K° = 1538 and r = 5x 10" 4 . 



Figure 25: (Q = (-8,8) 3 , u dQ = 0.004, y = y hex , fi = Aaiu) r*(f) for t = 
5, 10, 30, 50, 60; and T h (60) within O. Parameters are Nf = 512, N c = 32, 
K° = 98 and r = 10" 2 (left), and N f = 256, N c = 32, K° = 98 and r = 10" 1 
(right). 



ing errors. However, the same effect, on even the most refined 
meshes, can be achieved by adding random fluctuations to the 
model. In real life such fluctuations and changes in physical 
parameters are experienced by the growing snow crystal, as it 
moves through the atmosphere towards the earth. 

A numerical simulation with supers aturation uqq = 0.02 with 
J3 = /?fiat,3 is displayed in Figure|26] In this case capped columns 



appear, which can also be observed in nature; see I2ll35l. 

We end this subsection with computations, where we use a 
time dependent choice for ugQ. In particular, we set 



um(f) ' 



0.004 te [0,15) U [18,50], 
0.024 re [15, 18). 



(21) 



See Figure 27 for the results. First a plate forms and then, due 
to the fact that the supersaturation increases, the plate becomes 
unstable and new plate-like shapes grow at the corners of the 
plate. 

Finally, we perform two simulations, where we vary ft in 
time. In the first such example, we choose 



PmjiP) 'e[0,30), 
Pm,3(p) re [30, 50]. 



In a second example, we choose 



P{P) = 



£flat,l(^) 



t e [0, 20) , 
t e [20, 50] . 



(22) 



(23) 





Figure 27: (O = (-8, 8) 3 , u dQ as in (SlJ, y = y hex , (3 = /3 m , 3 ) T h (t) for t = 
15, 20, 30, 50; and 1^(50) within £1. Parameters are N f = 512, Af c = 32, 
^ = 98andr = 2xl0- 2 . 



Results for these choices of p and for uqq = 0.004 can be 



seen in Figure 28 We observe scrolls on plates, a shape that 
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Fi gure 28: (O = (-4,4) 3 , u dn = 0.004, y = y hex , J3 as in {22} (top), and as 
in (28) (bottom)) 1^(50). Parameters are N f = 128, N c = 16, ^ = 98 and 
t=10- ! . 

is also called plates with scrolls at ends, which also appear in 
the Magono-Lee classification of natural snow crystals 1361 , 
see also EH p. 46] and l35l . 

5. Conclusions 

We have demonstrated that an approach introduced by the 
authors in [QIOEEH provides a powerful computational tool to 
investigate pattern formation in crystal growth in a qualitative 
and quantitative way. The method makes it possible to simu- 
late facetted and dendritic growth simultaneously. We also ob- 
serve the instability of crystals leading to facet breaking. Many 
parameters in models for crystal growth are not known. The 
presented numerical method in combination with a compari- 
son to experiments can make it possible to estimate the relative 
sizes of parameters. In particular, by varying the condensation 
coefficient we were able to observe either plate-like growth or 
columnar growth. 

Let us finally summarise the results. 

• Surface energy effects taking anisotropy into account have 
been included in the model and, despite their small size, 
they turned out to totally change the character of the inter- 
facial dynamics. In our computations anisotropic surface 
energy is required in the interfacial dynamics to produce 
facetted dendritic growth. 

• The influence of the anisotropy in the condensation coeffi- 
cient, at least at small supersaturations, is not sufficient for 
facetted growth. 

• For small supersaturations the influence of the velocity 
term in ([7} is small in comparison with the curvature term. 

• The velocity at the tip of growing crystals depends in a 
linear way on the supersaturation. 

• Macroscopic models for crystal growth which are based 
on a diffusion equation in the gas phase, a mass balance 



on the vapour crystal interface and a modified Gibbs- 
Thomson law, taking attachment kinetics into account, are 
able to model a variety of phenomena in crystal growth, 
such as the appearance of solid plates, solid prisms, hollow 
columns, needles, dendrites, capped columns and scrolls 
on plates. 

Acknowledgement. The authors wish to express their deep 
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